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Abstract. The problem of the form of the 'arctic' curve of the six-vertex 
model with domain wall boundary conditions in its disordered regime is ad- 
dressed. It is well-known that in the scaling limit the model exhibits phase- 
separation, with regions of order and disorder sharply separated by a smooth 
curve, called the arctic curve. To find this curve, we study a multiple integral 
representation for the emptiness formation probability, a correlation function 
devised to detect spatial transition from order to disorder. We conjecture that 
the arctic curve, for arbitrary choice of the vertex weights, can be character- 
ized by the condition of condensation of almost all roots of the corresponding 
saddle-point equations at the same, known, value. In explicit calculations we 
restrict to the disordered regime for which we have been able to compute the 
scaling limit of certain generating function entering the saddle-point equa- 
tions. The arctic curve is obtained in parametric form and appears to be a 
non-algebraic curve in general; it turns into an algebraic one in the so-called 
root-of- unity cases. The arctic curve is also discussed in application to the 
limit shape of q-enumerated (with < q ^ 4) large alternating sign matrices. 
In particular, as q — > the limit shape tends to a nontrivial limiting curve, 
given by a relatively simple equation. 



1. Introduction 

In strongly correlated systems, the effect of boundary conditions can be relevant 
even in the thermodynamic limit. Consider for example a system whose parameters 
are tuned in such a way that at equilibrium it should be in a disordered phase, while 
its boundary conditions are chosen so that only ordered configurations are admis- 
sible in the proximity of the boundary. Due to the presence of strong correlations, 
it may happen that such boundary conditions induce ordered regions extending 
macroscopically from the boundaries deeply inside the bulk of the system. In such 
situation, spatial phase separation emerges, with ordered regions contiguous to the 
boundary, sharply separated from a central disordered region by a smooth curve, 
called arctic curve in the case of dimer models pQ . 

Essentially the same phenomena appear in other contexts, with different names, 
such as limit shape (in the statistics of Young diagrams |2j and rhombi tilings [3]), 
or interface (in random growth models for two-dimensional crystals [1]). More 
generally, the problem consists in finding the limit shapes and fluctuations of ran- 
dom two-dimensional surfaces arising, for instance, in plane partitions (or in three- 
dimensional Young diagrams, or in the melting of a faceted crystal) [SHZ], and also 
in dimer models on planar bipartite graphs with fixed boundary conditions, when 
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described in terms of the height function [HIE]. For recent developments see, e.g., 
[T0l[Tl] . In these contexts, the arctic curve is usually referred to as the frozen 
boundary of the limit shape. 

The standard example of an arctic curve is the famous arctic circle which ap- 
peared in the study of domino tilings of large Aztec diamonds [lj[12] . The name 
originates from the fact that in most configurations the dominoes are 'frozen' out- 
side the circle inscribed into the diamond, while the interior of the circle is a disor- 
dered, or 'temperate', zone. Further investigations of the domino tilings of Aztec 
diamonds, such as details of statistics near the circle, were also performed p~5HT5] ■ 

In the present paper we address the problem of the form of the arctic curve in 
the six- vertex model with domain wall boundary conditions. The case of generic 
Boltzmann weights of the disordered regime is considered. We find an explicit 
expression for the curve in a parametric form; in general, the curve appears to be 
non-algebraic. This property can be ascribed to the fact that the six-vertex model 
cannot be reduced to a model of discrete free (or Gaussian) fermions. 

Indeed, in all examples considered to date, the arctic curves appear to be alge- 
braic curves. At the same time, these examples can be seen as particular realizations 
of models of discrete free fermions, although on various types of lattice and with 
nontrivial boundary conditions. Some of them can be even reformulated as the 
six-vertex model at its free-fermion point with suitably chosen fixed boundary con- 
ditions (the correspondence being however usually not bijective). In particular, this 
is the case of domino tilings of Aztec diamonds, and the corresponding boundary 
conditions of the six- vertex model are exactly the domain wall ones [12j[16] . Thus 
the problem of the arctic curve extends naturally to the six-vertex model with 
generic weights, and with fixed (in particular, domain wall) boundary conditions. 
The phenomenon of phase separation for six-vertex model, in various regimes and 
with various fixed boundary conditions, was studied previously mainly numerically 
[T71 - H9] ; some analytical approaches to treat the problem were discussed in [20]]21"] . 

Historically, the six-vertex model with domain wall boundary conditions was 
introduced to prove the Gaudin hypothesis for norms of Bethe states [22J . The stan- 
dard framework for the model is the quantum inverse scattering method, invented 
in seminal paper [23j ; for a review and applications of the method see book [24] . 
The partition function of the model was obtained in terms of a determinant, known 
as Izergin-Korcpin formula [25 , 26] . The free energy per site was derived in [27]l28j. 
where phase separation was proposed as a possible explanation for the observed 
influence of boundary conditions on the thermodynamic properties. This, in turn, 
stimulated subsequent exact calculation of correlation functions of the domain-wall 
six- vertex model for generic values of its weights. This was done, using the quantum 
inverse scattering method, for one- and two-point boundary correlation functions 
|29|J30j , and, more recently, also for a particular non-local correlation function, the 
so-called emptiness formation probability [3T] . 

The interest in the domain-wall six-vertex model is also motivated by its close 
relationship with some problems in algebraic combinatorics. Apart from the already 
mentioned domino tilings of Aztec diamond, the model was also found to be related 
with enumerations of alternating sign matrices [12] . This observation was useful 
since it opened new possibilities in proving long-standing conjectures in this subject 
[32ll33j : see also book [34] for a review. In this context the arctic curve of the 
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domain-wall six-vertex model is of great interest since it describes the limit shape 
of very large alternating sign matrices [35l[36] . 

In the present paper, to address the problem of the arctic curve, we develop 
the idea proposed in [31] of studying the scaling limit of the emptiness formation 
probability, a correlation function with the capability of detecting spatial transition 
from order to disorder. Namely, in the scaling limit the emptiness formation prob- 
ability may have only a simple step-function behaviour, when varying coordinates 
from ordered to disordered regions, with the jump from one to zero occurring ex- 
actly at the arctic curve (actually, at one of its four portions, see Sections 2.3 and 
4.1 for details). This can be seen as a particular case of the general statement that 
the probability of finding a macroscopically large ordered sub-region must vanish in 
the disordered region. The arctic curve can then be obtained from certain multiple 
integral representation for the emptiness formation probability, as the condition on 
the relevant parameters (scaled coordinates) ensuring such a stepwise behaviour in 
the scaling limit. 

This programme was first fulfilled in [37] for a particular choice of the parame- 
ter A of the six-vertex model, the case A = 0, also known as the free-fermion point 
of the model. In that paper it is observed that the arctic curve can be obtained 
by investigating the multiple integral representation for emptiness formation prob- 
ability derived in |31j . It was found that the arctic curve corresponds to a rather 
peculiar solution of the system of saddle-point equations, namely, the solution char- 
acterized by a trivial Green function, with just one pole. In other words, the arctic 
curve appears to be in correspondence with the condition of 'condensation' of al- 
most all saddle points at the same value in the complex plane. In the case of 
A = this correspondence, between the condensation and the arctic curve, admits 
a rigorous derivation since in this case the saddle-point equations can be analysed 
using standard methods from the theory of random matrix models (although some 
complication arises, due to the two-cut nature of the problem, see Section 5.2, or 
also [37], for details). 

In the case of generic values of A, i.e., away from the free-fermion point, the 
system of coupled saddle-point equations describing the scaling limit of the empti- 
ness formation probability is extremely intricate and the standard tools of random 
matrix models (e.g., the Green function approach) are inapplicable. We conjecture 
that the arctic curve can nevertheless be found as the condition of condensation 
of almost all roots of the saddle-point equations to the same, known, value. As 
already observed in [36j . certain specific properties of the investigated multiple in- 
tegral representation, guaranteeing the step-function behaviour of the emptiness 
formation probability in the scaling limit, appear to be totally independent of the 
values of the parameters of the model, in particular, of the value of A. We also 
propose a method for deriving the arctic curve without a direct use of the Green 
function, but rather exploiting the fact that in the case of condensation the system 
of coupled saddle-point equations simplifies to a single equation (that we call 're- 
duced saddle-point equation'). This equation must necessarily have two coinciding 
real roots, for consistency with condensation itself. The arctic curve is then given 
as the condition of coincidence of two roots of the reduced saddle-point equation, 
where the value of this double root is the parameter which parameterizes the curve. 
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To obtain the explicit form of the reduced saddle-point equation, and hence 
to derive the arctic curve, one has also to find the thermodynamic limit of cer- 
tain function, which plays the role of generating function of a particular one-point 
boundary correlation function. Here we evaluate this limit for the whole disordered 
regime, thus extending the results of [36] where two relevant particular cases away 
from the free-fermion point were considered. We also give here a detailed exposition 
of the method. 

The organization of the present paper mainly follows the logic of our deriva- 
tion of the main result. After recalling the definition of the model, Izcrgin-Korepin 
formula, and the phase separation in Section 2, in Section 3 we introduce bound- 
ary correlation functions and address the problem of finding the contact points of 
the arctic curve with the boundaries, in the disordered regime. The definition of 
the emptiness formation probability, the multiple integral representation for this 
correlation function, and the corresponding scaling limit saddle-point equations are 
given in Section 4. We devote Section 5 to a derivation of the arctic curve at the 
free-fermion point of the model. Here we also explain how to derive the arctic 
curve without any use of the Green function but through the simpler condition 
of coincidence of two roots of the reduced saddle-point equation. In Section 6 we 
show that the two specific properties of the multiple integral representation for the 
emptiness formation probability which are relevant for the correspondence between 
the arctic curve and condensation are totally independent of the value of A. Here 
we also derive the reduced saddle-point equation and obtain the arctic curve in the 
disordered regime (|A| < 1); formulae (|6.17|) - (I6.19|) constitute the main result of 
the present paper. 

In Section 7, we discuss some particular cases of the main result, also in connec- 
tion with application of the model to the problem of ^-enumeration of alternating 
sign matrices (with < q ^ 4). In particular, we show that the limit shape 
of g-enumerated alternating sign matrices has a non-trivial limit as q — > 0; with 
a suitable choice of coordinates the corresponding limiting curve is just the cosine 
curve. In Section 8 we summarize the results and provide some concluding remarks. 

The text is followed by two appendices. In Appendix A we recall some results 
on the connection between the Izergin-Korcpin formula and the generating function 
of a particular one-point boundary correlation function. In Appendix B we obtain, 
for the disordered regime, the thermodynamic limit of this generating function, 
which is essentially used for the derivation of the main result. 

2. Six-vertex model, domain wall boundary conditions, and phase 

separation 

2.1. The model. The six-vertex model is canonically formulated in terms of 
arrows pointing along the edges of a square lattice. The arrows must obey the 
'ice- rule': there are two arrows pointing away from, and two arrows pointing into, 
each lattice vertex. 

Equivalently, one can describe the states of the edges in terms of variables 
taking two values, e.g., 0, 1. We shall use the convention that such a variable takes 
value (1) if the arrow on the edge is pointing upward or right (downward or 
left). Let us consider a single vertex of the square lattice and let fi, v, v 1 and fi' 
be the variables attached to the left, bottom, top, and right edges of the vertex 
respectively. Let w(fi, v, z/, fi') denote the Boltzmann weight associated with the 
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vertex. In the six-vertex model 

w(fi, v, z/, fjf) = if ii + v^n' + v'. (2.1) 
Imposing the arrow-reversal invariance of the Boltzmann weights, we set 

io(0,0,0,0) = io(l,l,l,l) = a, 

w(0, 1, 1, 0) = io(l, 0, 0, 1) = 6, (2.2) 
u;(0,l,0,l) =io(l,0,l,0) =c, 

where a, b, and c are some real positive constants, referred often below as 'weights'. 
The partition function of the six-vertex model can be defined as 

Z = II /'-/•"'...•"'.,■ •/'!,,*• ( 2 - 3 ) 

Here the product is taken over vertices of the lattice and the sum is performed over 
values of edge variables; the variables are subject to the linking conditions 

= /'< I ..- v'ij = (2.4) 

To define completely the model, one has also to specify the boundary conditions 
in (|2.3|) . We shall consider the so-called domain wall boundary conditions, defined 
below. 

From the analysis of thermodynamic limit of the model with periodic boundary 
conditions (see, e.g. [38] ) it is well known that the model has three physical regimes, 
or phases; two of them are regimes of order and one is a regime of disorder. Let us 
introduce the parameter 

* = t- 

The case A > 1 is called ferroelectric regime. Depending on whether a > b or 
a < b, ground state configurations are built from vertices of weight a, or b; in 
these configurations all arrows are ordered ferroelectrically. The case A < — 1 is 
called anti-ferroelectric regime. In this case, the doubly-degenerated ground state 
is formed only by vertices of weight c, and thus the directions of arrows alternate 
along vertical and horizontal rows of the lattice. Finally, the case | A| < 1 is called 
disordered (or critical) regime. In this case typical configurations contain vertices 
of all three weights; there is no particular order of arrows in these configurations. 
All over the paper we use the following parameterization for the weights: 

a = sin(A + 77), b = sin(A — 77), c = sin 2r). (2.6) 

We also have 

A = cos 277. (2.7) 

The parameter A has meaning of a rapidity variable and 77 is the so-called crossing 
parameter. 

In the case of the disordered regime, i.e., when a, b, and c are such that | A| < 1, 
both A and 77 are real and satisfy 

tKA<7t-77, 0<?7<|. (2.8) 

The other two regimes, the ferroelectric and anti-ferroelectric ones, can be ap- 
proached by analytic continuation (modulo a purely imaginary common factor for 
each weight), choosing A = iA, A = iff and A = n/2 + iA, A = n/2 + ifj, respectively, 
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where A and fj are real, and satisfy some further restrictions to ensure positivity of 
the weights. 

For later reference, let us briefly discuss some symmetries of the model. First, 
we mention the so-called crossing symmetry. It is a symmetry of the Boltzmann 
weight under reflection of the vertex with respect to the vertical (horizontal) line, 
reversing the values of the edge variables on the horizontal (vertical) edges, and 
exchanging the values a <-> b. For example, the reflection with respect to the 
vertical line gives the crossing symmetry relation 



w(n, v, v' , fi'\\) = w{jjl , v, v' , [i\-k — A). 



(2.9) 



Here the bar means that the value of the edge variable must be reversed, i.e., p, — 1 
(p, = 0) if /x = (// = 1). Similarly, one can also consider the reflection with respect 
to the horizontal line, that gives an analogous relation. 

Another useful, and in fact even simpler, property of the six-vertex model 
Boltzmann weight is its invariance under reflection of the vertex with respect to 
one or another diagonal, 



w(p, v, v , /i'|A) = w(v, /i, //, v'\ A) = wiy' , //, /i, v\ A). 
We shall refer to this property as the diagonal-reflection symmetry. 



(2.10) 



2.2. Domain wall boundary conditions. Consider the six-vertex model on 
a square lattice obtained by intersection of an equal number of vertical and hori- 
zontal lines. Domain wall boundary conditions correspond to hxing the variables 
on the external edges of the lattice in a specific way, namely, with all arrows on 
horizontal edges pointing outward the lattice, and all arrows on external vertical 
edges pointing inward the lattice. Equivalently, these boundary conditions mean 
that all variables attached to external edges on the top and on the left are set to 1 , 
while those attached to external edges on the bottom and on the right are set to 0. 

Denoting by N the number of vertical, or, equivalently, horizontal lines of the 
lattice on which the model with domain wall boundary conditions is considered, we 
denote by Zn the partition function of this model, 



N 

Zn = II W ^ ' Vi >3 ' V 'i,3 ' 



(2.11) 



Here the stars indicate that the subscripts must run over all possible values, e.g., 
Hn,* = 1 means that (j,n,i = MJV.2 • • • = hn,n = 1. 

In parametrization (|2.6p . it can be shown that the partition function admits 
the following exact representation [25] 



1 jV 



[ sm(A - 77) sm(A + 77)] 
Here Dn = Dn(A) stands for the Hankel determinant 



(2.12) 



D 



N 



JN-l) 



(A) 
^W(A) 



p(*-D(A) ^W(A) 



.(2JV-2) 



(A) 



(2.13) 
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where function ip(X) is given by 

sin 2?y 

^ (A) := sin(A-r / )sin(A + ?? ) - (2 ' W) 

Formula (|2.12[) follows from the more general result, known as Izcrgin-Korepin 
formula (see Appendix A, equation l|A.2p ). which is valid for the model with inho- 
mogeneous weights |25] . The proof of Izergin-Korepin formula is essentially based 
on the quantum inverse scattering method [521[Mj; see also for details. 

Using formula (|2.12[) one can study the partition function in the thermodynamic 
limit. The quantity of interest is the free energy per site /, defined as 

/ := _ um (2.15) 

It can be shown [2"TU2"8"ll3T)] that the free energy per site of the model in the disordered 
regime is given by the expression 

/ = _ ln / QS in(A + ^)sin(A-r,)\ 
\ sina(A — rj) ) 



where we have used the notation 



(2.17) 



7T — 2?7 

For a derivation of formula (|2.16p see also Appendix B. 

Formula (|2.16[) was obtained for the first time in [27], where also the case of 
ferroelectric regime was considered; the case of the anti-ferroelectric regime was 
studied in [28 . In these studies, the following observation have been made: the 
free energy per site is greater in the case of the domain wall boundary conditions, in 
comparison with the case of periodic boundary conditions, both for the disordered 
and anti-ferroelectric regimes, but remains the same for the ferroelectric regime, 

/domain wall > /periodic, for A<1. (2-18) 

This result can be explained by the fact that in these two regimes, admissible con- 
figurations are significantly constrained by the domain wall boundary conditions. 
Moreover, for large lattices, these effects are macroscopically large, i.e., the model 
exhibits phase-separation phenomena. 

2.3. Phase separation and arctic curves. The six-vertex model exhibits 
spatial separation of phases for a wide choice of fixed boundary conditions [17j . 
Roughly speaking, the effect is related to the fact that ordered configurations on the 
boundary can induce, through the ice-rule, a macroscopic order inside the lattice. 
Analytically, as already mentioned, a signal for presence of phase separation is a 
change of the free energy per site, in comparison with the case of periodic boundary 
conditions. 

The notion of phase separation acquires a precise meaning in the scaling limit. 
This is a thermodynamic limit which is to be treated as a continuum limit, namely, 
in this limit the number of lines of the lattice (in each direction) tends to infinity 
and the lattice spacing vanishes, while the total size of the system (in each direction) 
is kept fixed. 

For the domain- wall six- vertex model we take N — > oo, but the whole lattice 
is scaled to a finite square, e.g., with sides of length 1. Typical configurations are 
constrained to have macroscopic regions of ferroelectric order near the boundary. 
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Figure 1 . The four regions of ferroelectric order and the region 
of disorder for the model in the scaling limit, in the disordered 
regime. 

More precisely, if the parameters of the model are tuned to the ferroelectric regime 
then there is full compatibility between boundary and bulk phase, no conflict arises, 
and the whole system is in the ferroelectric order. If, however, the parameters 
of the model are tuned to the disordered or the anti-ferroelectric regimes, then 
there is a competition with the ferroelectric ordering induced by the boundaries, 
and phase separation emerges. The present understanding of the phase-separation 
phenomenon in these two regimes is mostly due to numerics |18U19| ; theoretical 
considerations based on a variational principle approach can be found in [20H21j . 
Below we describe the case of disordered regime, and mention briefly the situation 
in the anti-ferroelectric regime at the end. 

In the case of the disordered regime, in the scaling limit five regions appear: 
four regions of ferroelectric order, Fnw, Fne, Fse, Fsw in the four corners of the 
square, and one region of disorder D, in the centre, see figure [TJ The region of 
disorder is sharply separated by some smooth curve A, called the arctic curve. The 
arctic curve and the square have four contact points, located each one on a side of 
the square. 

The arctic curve in the disordered regime consists of four portions, 



where i -1 * separates the region (i — NW, NE, SE, SW) from the internal region of 
disorder D. 

To describe, for example, the curve -Tnw, one can introduce a function T(x, y; A) 
where x and y are coordinates in the scaling limit and A is the parameter of the 
weights, 



It is to be emphasized, that the explicit form of T(cc, y; A) is significantly determined 
by the value of r\. Given function T(x, y; A) one can readily obtain all other portions 



A — 2~nw u r NE u r'sE u -Tsw, 



(2.19) 



r N w : A) = 0, 



x,y £ [0, k}. 



(2.20) 
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of the arctic curve. Due to the crossing symmetry, relation (|2.9|) . we have 

r NE : T(l-x,y;n-X) = 0, x G [1 - K, 1], y G [0, k]. (2.21) 

Exploiting the diagonal-reflection symmetry, we further have 

r SE : T(l-z,l-2/;A) =0, x,y G [1 - k, 1]. (2.22) 

Finally, and essentially similarly, we have 

r sw : T(z,l-y;7r- A) = 0, a: G [0, k], y G [1 - k, 1]. (2.23) 

Thus, knowledge of function T(x,y]X) allows one to determine the whole arctic 
curve A. 

It turns out that function T(x,y;A) is algebraic (for example, quadratic in 
x and y) only for special values of rj, while in general it is a rather complicated 
non-algebraic, or transcendental, function. To describe the curve -Tnwj m practical 
calculation, we shall use the parametric form 

x = x(0, y = Y((), Ce[0,Co], (2.24) 

where X((), Y(Q and Co depend on the parameters of the model, A and r\. This 
must be taken into account when reconstructing the whole arctic curve A as ex- 
plained above. Functions X(Q and Y(() will be monotonously decreasing and 
increasing functions, respectively, satisfying X(0) = k, Y(0) = and X(fo) = 0, 
Y"(Co) = K- Furthermore, these two functions are in fact just a single function, since, 
due to reflection symmetry, we have the relation X(£) = Y((q — £). The explicit 
results are given in Section 6, and some interesting particular cases are discussed 
in Section 7. 

To conclude here, we mention briefly that in the case of the anti-ferroelectric 
regime the picture of phase separation, though reminding the situation of the disor- 
dered regime, is more complicated. The competition between the anti-ferroelectric 
central region and the ferroelectric boundary regions induces the emergence of an 
intermediate disordered region. As a result, there are two phase-separation curves: 
an outer one, which is the usual arctic curve, separating the ferroelectric and dis- 
ordered phases, and an inner one separating the disordered and anti- ferroelectric 
phases. We refer for further details to [THIUS] where numerical simulations for this 
regime were performed. 

3. Boundary correlation functions 

3.1. Definitions and properties. Before addressing the whole problem of 
the arctic curve, it is useful to restrict first to the simpler problem of the position of 
the contact points. This problem can be treated by studying correlation functions 
describing probabilities of certain configurations near the boundaries, the so-called 
boundary correlation functions. 

Following |29j , we consider two kinds of boundary correlation functions, denoted 
as and G$ (in these notations the superscript in parenthesis, r, refers to 

some distance on the lattice, and should not be confused with an rth derivative 
with respect to any variable). As it will become clear later, the first correlation 
function will play a somewhat fundamental role, through its generating function 
/ijv(^) which enters various correlation functions, while the second one is more 
specific but relevant for addressing the problem of the contact points. 
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(r) 

The first correlation function, , is designed to reflect the fact that the con- 
figurations of the model with domain wall boundary conditions admit one and only 
one c-weight vertex on each of the four boundary lines of the lattice. Specifically, 

(r) 

we define H)^ as the probability that the sole c-weight vertex of the first horizontal 
line from the top occurs at the rth position from the right. It can be introduced 
by the formula 



N 

where S(p, //) stands for the Kronecker symbol, 



(3.1) 



[0 if fi^ n' 
This correlation function, satisfies the sum-rule 

N 

= (3.3) 

r=l 

which simply expresses the fact that the total probability of finding a c-weight 
vertex anywhere in the first horizontal is exactly one. 

Some properties of correlation function can be obtained using the sym- 
metries of the model. In particular, crossing symmetry implies 



<(A)=<-' r+1 V-A). (3-4) 

(r) 

Using the diagonal-reflection symmetry, one finds that function H y N ' also gives, for 
example, the probability that the sole c-weight vertex on the rightmost vertical line 
occurs at rth position from the top. 

(r) 

The second correlation function, G N , is defined as the probability of finding 
a given state on a horizontal edge of the first line, i.e., it is a particular case of 
polarization. For r = 1, . . . , TV — lwe define this correlation as the probability of 
having the edge variable fixed to one, on the edge connecting rth and (r + l)th 
vertices (from the right) of the first horizontal line (from the top), 

N 

G N = Z n X <KMr,i,l) Yl ^(MiJ.^i.^i.Mij) 

M,M i,j=i 



, ■ (3-5) 

^*,N=A I i »=0 



We can also set G { x ] = 1 and G { ° } = since the edge variables on the leftmost 
and the rightmost external edges are already fixed by the boundary conditions to 
fJ-N,i — 1 and n'n = 0. 

As noticed in [29], correlation functions H N and G N are in fact closely related. 
Indeed, since 5(/x r ,i, 1) + 5(/x ri i,0) = 1, formulae (|3.ip and f|3.5|) imply that Gpj = 
H n ] + G iv _1) - Taking into account that G^ } = 0, one can equivalently bring this 
relation into the form 

G$=Htf + H%-V+... + H$\ (3.6) 

Recalling relation f|3. 3|) . we obtain that (|3.6p implies that G^ = 1, in agreement 
with the definition of this correlation function. 
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(r) 

As we shall see below, in the scaling limit function G N turns into a step 
function; the point where its value jumps from to 1 is exactly the contact point 
of the arctic curve. Relation (|3.6|) ar 
us to find the location of this point. 



(r) 

of the arctic curve. Relation (|3.6[) and certain properties of function H N will allow 



3.2. Function h^(z) and its large N limit. The correlation function 
can be cquivalcntly considered through its generating function 

N 

M*) := X^&V" 1 - ( 3 - 7 ) 

r=l 

Due to (|3.3p . the generating function satisfies the normalization condition /ijv(l) = 
1. For later referring, let us also mention the property 

h N (0)=HV = Ca2N " Z *-\ (3.8) 

This property simply follows from the fact that the c-weight vertex appearing in 
the corner induces domain wall boundary conditions on the remaining (N — l)-by- 
(N — 1) sub-lattice [40 . We have also the relation 

h N {z-\)=z N - 1 h N {z- 1 -ii-\) 1 (3.9) 



which is a direct consequence of relation (|3.4|) . 

Given function tiN(z), and taking into account relation (|3.6j) . one can readily 

(r) 

represent correlation function G N as follows 



<# = -^* 7^7 dz. (3.10) 



,0) _ J_ f Mf) 

2vri J Co (z - l)2f 

Here Co is a simple closed counterclockwise-oriented contour surrounding point 
z = 0, and lying in its small vicinity. In the next Section, a much less trivial 
example of correlation function will be given, where a multi-variable generalization 
of function hisr(z) will appear. 

In |29j . where the boundary correlation functions were introduced, it was shown 
that they can be expressed similarly to formula (|2.12[) for the partition function, 
but with properly modified entries in the last column of the determinant. Here 
we rely on a similar but actually different representation, which relates generating 
function h^(z) with the partition function of the 'partially' inhomogeneous model 
with a single parameter of inhomogeneity. This representation is discussed in detail 
in Appendix A, see equations (|A.8|) and (|A.12|) . 

The following result plays a crucial role for what follows. 

In the disordered regime, i.e., with A and r/ satisfying (|2.8|) , the function 
ln/ijv(z), for z real and non- negative, has the following N — > oo behaviour 

1*M7(0) =Nln( ^-^^-^f.^ ) + o(N), (3.11) 
\ a sm( A — rj) sin a(f + A — r]) sin f / 



where 



(n _ sin(A + rj) sin(g + A - rj) 
7lU " sin(A-77)sin(£ + A + T7) [ ' 



and a has been defined in (|2.17D . 

The proof of (|3.1ip is given in Appendix B. 
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3.3. Contact points. To find the contact points in the disordered regime, 

(r) 

i.e., the value of constant k, we study correlation function G N in the scaling limit. 
Let us define function G(x) (0 ^ G(x) ^ 1) as the scaling limit of correlation 
function Grfc, 

n N — r 

G(x):= lim G%\ xe[0,l]. (3.13) 

r,N—>oo 1\ 

In the scaling limit, the dominating configurations are such that all edge variables 
for horizontal edges of the first line are fixed to 1 in the region Fnwj an d to in 
the region Fne- Therefore, G(x) must have a simple stepwise behaviour, with the 
jump occurring precisely at the contact point, 

G(x) = l 1 f ° r °^< K . (3.14) 
[0 for k<^1 

Hence, the value of k can be found by addressing the question of how correlation 

(r) 

function G N may exhibit a stepwise behaviour, as N, r — ► oo. 

To address this question, we consider integral representation (|3.10p , and apply 
the saddle-point technique to analyze its behaviour in the scaling limit. Formally, 
the saddle-point equation reads 

AT^oo N J 

In principle, to investigate the solution of the saddle-point equation, one would need 
to know the last term in (|3 . 1 5[1 for z everywhere in the complex plane. However, 
it turns out that in addressing the problem of contact points the complete saddle- 
point analysis can be avoided, once the mechanism at the origin of the stepwise 
behaviour of function G(x) is understood. Furthermore, formula (|3.11| . which was 
derived under the assumption of z being real and positive, appears to be sufficient. 

Indeed, let us first restrict to the case A = 0, where function h^(z) is given by 
asimple explicit expression, h^{z) — [|(1 — sm2A)2+|(l + sin.2A)] JV ~ 1 , see [291141] 
for a proof. From this expression it follows that the saddle-point equation has only 
one solution, which we denote by z s p It can be seen that z s p which depends on 
the value of x £ [0, 1], is always real and positive, with the corresponding steepest- 
descent contour perpendicular to the real axis. When deforming the contour Cq in 
(|3.10P to the steepest-descent contour, the position of the saddle-point z 3 . p . with 
respect to the pole at z = 1 of the integrand appears to be crucial. Namely, for 
z s .p. < 1, the pole at z = 1 can simply be ignored in the deformation of the contour, 
and the saddle-point evaluation of the integral vanishes in the scaling limit. On 
the other hand, for z s p . > 1, in deforming the contour through the saddle-point, 
we necessarily cross the pole at z = 1, thus picking a contribution equal to minus 
the residue of the integrand at this pole. This mechanism is at the origin of the 
stepwise behaviour of G(x), with the step occurring when x is such that z s . p . = 1. 
This last condition thus determines the position of the contact point. 

When considering the problem for generic values of A, we note that hjy(z), 
being a polynomial, is analytic, and the only pole to be taken into account in the 
integrand of (|3.10[) is again the one at z = 1, with residue —1, due to the property 
/ijv(l) = 1- Moreover, it can be verified that saddle-point equation (|3. 15|) has only 
one solution when z is real and positive, as assumed in formula (|3.11| . The stepwise 
behaviour of G(x) can thus be ascribed again to the fact that, for some particular 
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value x = k € [0,1], this real and positive saddle-point occurs at z = 1. This 
immediately determines the value of the contact point, 



i , ,. lnh N (z) 
K = 1 - | lim — 

N^oo N 



(3.16) 



It is apparent that, for the evaluation of the contact point, the knowledge of the 
large N behaviour of ln/ijv(z) for z real and positive, formula (|3.1ip . is sufficient. 

Prom expression (I3.16|) . and asymptotic behaviour (|3.11|) . in the case of generic 
value of A in the disorder regime, we obtain the following explicit expression 



a cot a(X — i]) — cot(A + 77) 



(3.17) 



cot(A — rf) — cot(A + 77) 

Note that, under the replacement A t— > 7r — A, we have k i— > 1 — k, as expected, due 
to the crossing symmetry. In particular, for X — n/2 we have k = 1/2 for all values 
of rj in the disordered regime. 

4. Emptiness formation probability 

4.1. Definition and basic properties. We want now to address the problem 
of finding the arctic curve A for the domain-wall six- vertex model. The most direct 
way would be to consider a one-point correlation function, such as polarization of all 
horizontal (or vertical) edges on the lattice, compute for it a suitable representation, 
and finally investigate its behaviour in the scaling limit. This would give access not 
only to the arctic curve A but to the whole 'limit shape' of the model, in its 
'height function' formulation, in the sense of [21]. The main obstruction for such 
strategy resides in the difficulty of computing explicitly such one-point correlation 
function. If one restricts to the arctic curve, however, it is possible to devise 
some other correlation function, maybe rougher than polarization, but nevertheless 
refined enough to discriminate the spatial transition between order and disorder, 
and simpler to calculate. 

Indeed, consider the emptiness formation probability (EFP), namely, the prob- 
ability F^' sS> of having the edge variables fixed to 1 on the first s horizontal edges, 
from the top of the lattice, and located between the rth and (r + l)th vertical lines, 
from the right [3T] 

f n' S) = z n 1 X! < 5 (^r.i,l)<5(AV,2,l)-"<5(AV, ;5 ,l) 

N 



, ■ (4.1) 

i / »,w=Mi I *=0 

Evidently, EFP is just a particular case of a general s-point correlation function, 
but it enjoys specific properties, outlined below, which are suitable to address the 
problem of the arctic curve of the model in the disordered regime. 

By construction, EFP as defined in (|4.ip actually measures the probability that 
the edge variables of all edges in the top-left (N — r) x s sub-lattice are fixed to 
1. This follows directly from the definition of the six- vertex model, see condition 
(|2.1|) . and the domain wall boundary conditions. In other words, EFP measures 
the probability that the top-left (N — r) x s sub-lattice has ferroelectric order, 



14 



F. COLOMO AND A.G. PRONKO 



To find the arctic curve, one must address the asymptotic behaviour of EFP in 
the scaling limit. Namely, we are interested in the limit N, r, s — > oo, while keeping 
the ratios fixed. We set 

N — r s 
x:=^ r , y-=^> y e [0,1]. (4.2) 

In this limit, coordinates x and y will parameterize the unit square to which the 
lattice is rescaled. Correspondingly EFP is expected to approach some limiting 
function 

F(x,y):= lim Fp s) . (4.3) 

N ,r,s — >oo 

Since EFP measures the probability that the top-left (N — r) x s sub-lattice 
has ferroelectric order, in the scaling limit, F(x, y) should tend to 1 whenever 
(x,y) G Fnw- On the other hand, in the disorder region, by definition, the number 
of edges breaking the ferroelectric order is macroscopic. Thus, as soon as (x, y) 
enters the disordered region D, function F(x, y) should immediately vanish. As a 
result, in the disordered regime, the limiting function is just the step function, 

F(xy) = { 1 if ^ gFnw (44) 
[0 if x,ye D U F NE U F SE U F sw ' 

where we have used notations of figure [1] 

Hence, EFP can be used to find curve Tnw, which separates region Fnw from 
the remaining part of the unit square. We recall, see Section 2.3, that knowledge of 
curve -Tnw allows one to recover other portions -Tne, -Tse, and Isw, thus solving 
the problem of the arctic curve of the model. 

Curve -Tnw can be found by analysing a suitable representation for EFP in the 
scaling limit. Moreover, formula (|4.4p implies that in such an analysis, similarly to 
the case of the contact point considered in Section 3.3, understanding the mech- 
anism at the origin of the stepwise behaviour of EFP in the limit is sufficient for 
finding the curve. 

4.2. Multiple integral representations. In the following we resort to the 
representation for the emptiness formulation probability in terms of s-fold integral 
obtained in [31]. The calculation itself is based on the quantum inverse scattering 
method [23, 24J and it represents further development of ideas of [29l[42] where 
boundary correlation functions were studied. In fact, for what follows we need two 
equivalent multiple integral representations, which differ only in their integrand, 
one being just the totally symmetrized version of the other. 

To give these representations we need to introduce some notations first. For 
s = 1, . . . ,N, let us define functions hi\r lS (zi, . . . , z s ), where the second subscript 
refers to the number of arguments, by the formula 



in,s(zi, ■ ■ ■ , z s ) := TT 



Zk 
l<i<fc<s 



z{ 1 h N - a+ i(zi) ... z s s 1 h N _ s+1 (z s ) 

~ 2 (zi - l)h N - s+2 (zi) ... z s s ~ 2 (z s - l)h N - s+2 (z s ) 



(4.5) 
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Note that we also have fiN,i(z) = hxlz). The functions flN,s(zi, ■ ■ ■ ,z s ) are sym- 
metric polynomials of degree (N — 1) in each variable z\, . . . , z s . In fact, these 
functions can be recursively constructed starting with function fiN,N(zi, ■ ■ ■ ,Zn), 
since, due to the structure of (|4.5p . it is easy to observe the relation 



h NtS (zi, . .., z s -i, 1) = h NtS -i(zi, . . . , « a _i). (4.6) 

One more relation of the same kind is 

h N>s (zi, z„-i, 0) = h N (0)hff-i, a -i(zi, . . . , z a -i). (4.7) 

where the value of /ijv(0) is given by formula (|3.8|) . As shown in [31^143] . functions 
hN,s{z\, . . . , z s ) are closely related to the partition function of the partially inho- 
mogeneous six- vertex model with domain wall boundary conditions, with s nonzero 
inhomogeneities. This relationship is sketched in Appendix A. 

We are now ready to turn to the multiple integral representations. The first 
multiple integral representation reads: 

(r>s) _ {-IT r f yr [(t*-2At)z j + l]°-* 

x n t ^z k z -~2Atz J+ i hN -^-- zs)dzi --- dzs - (4 - 8) 

Here Co, as before, denotes a simple anticlockwise oriented contour surrounding 
the point z = and no other singularity of the integrand. We use, besides the 
anisotropy parameter A, the asymmetry parameter 

t := -. (4.9) 
a 

The second, essentially equivalent, representation reads: 

F (r,s)_ j-iy^wz. r r A [ft 2 -2Af)z J + ir 1 
x TT -5 K TT (zk-zif 

J- 1 t 2 z jZk - 2AtZj + 1 11 v 3 ' 

j.k=l 3 K 3 l^j<k^s 

x h,N,s{zi, ■ • ■ , z s )h StS (ui, . . . , u s ) dzi ■ ■ ■ dz s . (4.10) 

Here Z s denotes the partition function of the six-vertex model with domain wall 
boundary conditions on an s x s lattice, and 

Uj := — — Zj ^- . (4.11) 

1 (t 2 - 2At)zj + 1 y ' 

The integrand of (|4.10p is just the symmetrized version of the integrand of (|4.8[) . 
with respect to permutations of the integration variables z\,...,z s - This follows 
through certain symmetrization procedure, see Appendix D of [44] . and an addi- 
tional identity proven in [31] . 

Representations (|4.8|) and (|4. 1 0(1 hold for all the three regimes of the model. 
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4.3. Saddle-point equations. We are interested in the behaviour of EFP in 
the so-called scaling limit, that is in the limit where r, s and N are all large, with 
the ratios r/N and s/N kept finite (and smaller than 1). In this limit, in principle, 
EFP can be analysed through the saddle-point method applied to multiple integral 
representation (|4.10l) . 

To find the corresponding system of saddle-point equations, we apply standard 
arguments (familiar, for example, from the matrix model context) to the integrand 
of (|4.10[) . We have to use the fact that quantities like ln/ijv,s(zi, ■ ■ ■ , z s) are of order 
s 2 , and that their derivatives with respect to Zj's are of order s. As a result, we 
arrive at the following system of coupled saddle-point equations 

s r s(t 2 - 2At) -A / t 2 z k - 2At 
Zj - 1 Zj {t 2 - 2At)zj + 1 ^ \ t 2 ZjZ k - 2AtZj + 1 

Pzk 2 \ din h NtS (zi, ...,z s ) 

t 2 ZjZk — 2Atzk + 1 Zk — zj J dzj 

t 2 - 2At + l d\nh s ^( Ul ,...,u s ) = 
[{t 2 -2At)z 3 + l] 2 duj 

where j = 1, . . . , s. In writing this expression we neglect all sub-leading contribu- 
tions (i.e., estimated as o(s)). 

As we shall see, although the saddle-point analysis of the multiple integral 
representation for EFP cannot be actually performed in the full rigour in general, 
nevertheless we will be able to find the condition on the ratios r/N and s/N which 
corresponds to the jump from to 1 of EFP in the scaling limit, i.e., to derive the 
equation for the arctic curve. In many respects, the whole procedure will remind 
the saddle-point analysis of the single integral (|3 . 1 0[) which allowed us to find the 
location of the contact point. 



5. An example: the case A = 

5.1. Preliminaries. Before to proceed with the general case, it is useful to 
consider technically the simplest case A = 0. The approach described in this 
section was given in detail in [37] . Our aim here is to recall the main steps of the 
approach and formulate a modified version of it, which appears sufficiently simple 
and powerful to be applied to the much more complicated case of generic values of 
A in the disordered regime, |A| < 1. 

In the case A = 0, or, equivalently, rj — 7r/4, function hpf(z) is known explicitly. 
It has a very simple form (see, e.g., [4T] ) 

/t 2 z + l\ N - 1 

kN( - z)= {^TT) ' i = tan(A-f). (5.1) 

The determinant entering the definition of functions hi\r, s {zi, . . . , z s ) turns out to 
reduce simply to a Vandermonde determinant; straightforward calculation gives us 

j=l V ' l^j<k^s 
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From this formula, using Uj = (1 — Zj)/(t Zj + 1), we also have 

s 1 

h a , a (u u ...,u B ) = J[ (t 2 z , 1 vs-i II ( t2z j z k + l)- (5.3) 

Taking into account that in the considered case Z s — 1 (see, e.g., [55]) and hence 
Z s /a s ^ s ~ x ^c s = (t 2 + 1) S ( S - 1 )/ 2 J we obtain, as a result, that representation (|4.10j) at 
A = reads: 

(r , s) (-I)*-*)/* r r * (t* Zj + !)»-» 



x JJ (zj - z fe ) 2 dzi • • • dz s . (5.4) 

Due to the last factor of the integrand, which is the squared Vandermonde determi- 
nant, formula (|5.4p naturally recalls the random matrix model partition function, 
that will be exploited in what follows. 

Here we would like also to mention that the last expression can be brought into 
the form 

F^' s) = det Q (5.5) 
where matrix Q is an s-by-s matrix with the entries 

Here is an arbitrary parameter whose value does not affect the determinant in 
(|5.5|) . Note that matrix Q has the structure of a Toeplitz matrix, i.e., its entries 
depend on the indices only through their difference. 

The determinant formula (|5-5[) has the following important implication. Let us 
consider the quantity 



r(r,s) _ (-i)'t^)/ a r r A(^+iT 

N ■ s !(27ri)«(^ + l)^- S )/ c -''7 cr il z^Zj-iy 



x J| (2 J --2 fe ) 2 dz 1 ---(k s . (5.7) 

This quantity differs from (|5.4p only in the integration contours: here C-f is a closed 
clockwise oriented contour (the minus sign in the superscript indicates negative 
direction) in the complex plane enclosing point z — 1, and no other singularity of 
the integrand. Since the integrand is kept intact, we have 



= detS, (5.8) 



where entries of matrix S differ from those of matrix Q only in the integration 
contour, 

Sjk = , 9 1 nv i 1)N ~ S (z - ey-^- 1 dz. (5.9) 

3 2iaif + 1) N - S Jc- z r {z-l) s v ' y ' 

Here is again an arbitrary parameter, whose value does not affect the determinant 
in (|5.8[) . Setting 9=1, matrix S reduces to an upper-triangular matrix, with all 
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its diagonal entries equal to 1. Hence 

Ip s) = 1, (5.10) 

for all r, s = 1, . . . , N. Identity (|5.10[) . together with some other properties of repre- 
sentation (|5.4[) . turns out to be one of the main ingredients entering the derivation 
of the arctic curve. 

5.2. Matrix model approach. We want now to address the asymptotic be- 
haviour in the scaling limit, as defined in Section 4.1, of multiple integral repre- 
sentation (|5.4[) for EFP in the case A = 0. Inspired by the obvious analogy with 
the large s limit of the partition function of an s x s Random Matrix Models, we 
rewrite multiple integral representation ()5.4j) as follows: 

Fir , s) _ (-ir (s+1)/2 / 

N S !(l+i 2 )(i/y-i)s 2 (27ri)s J Co ' 




z k ydz 1 ---dz s (5.11) 



where the potential is given by 

V(zj) = ln(Zj - 1) + - \nzj- ( - - 1 J In^zj + 1). (5.12) 
y \ y J 

Here we have used the scaling variables x and y introduced in (|4.2p . The corre- 
sponding system of coupled saddle-point equations reads 

1 1 - x (1 - y)t 2 2 -A 1 



1 y Zj y{t 2 Zj + 1) s^zj- z k 



-E^— = ' ( 5 - 13 ) 

k=\ 



where j = 1, . . . , s. 

The standard physical picture reinterprets the saddle-point equations as the 
equilibrium conditions for the position of s particles, each with charge 1/s, with 
logarithmic electrostatic repulsion, in an external potential. In the present case the 
latter is built from three logarithms, and can be seen as generated by three external 
charges, 1, (1 — x)/y, and — (l/y — 1) at positions 1, 0, and —1/t 2 , respectively. It 
is natural to refer to this model as the triple Penner model. Although the simple 
Penner [45j matrix model has been widely investigated, not so much is known about 
the more complicate multiple or generalized Penner models. The main ideas we 
shall use in the following are based on the phenomenon of eigenvalues condensation, 
typical of Penner models, and nicely described in [4T>ll4"T] . 

To investigate the structure of solutions of the saddle-point equations (|5.13|) 
for large s, one can start with introducing the Green function 

£?(*):= i£-J_, (5.14) 

s z — Zj 

3 = 1 J 

which, if the Zj's solves (I5.13|) . has to satisfy some particular differential equation, 
which can be derived by standard means. In the scaling limit, such differential 
equation reduces to an algebraic one for the limiting Green function. In principle, 
the latter can be solved, and an expression for Qiz) can be obtained in the large 
s limit. In this limit, the poles of Q(z), organize into cuts in the complex plane, 
and the discontinuities across such cuts define, when real and positive, the density 
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of solutions of the saddle-point equation as s — * oo. In the present case, however, 
the Green function Q{z) is not completely determined by the above procedure, 
since it contains an extra parameter, which can be fixed only implicitly. This is 
a direct consequence of the 'two-cut' nature of the problem (see [37] for details). 
The problem of explicitly finding the density of solutions of saddle-point equation 
(|5.13| as s — > oo, for generic values of x, y, is therefore a formidable one, not to 
mention the evaluation of the corresponding 'free energy', and of the saddle-point 
contribution to the integral in (|5.11|) . 

The problem we are addressing is fortunately much more modest: we are 
presently interested only in the expression of curve -Tnw, which, as explained in 
Section 4.1, corresponds to the curve in the unit square, where, in the scaling limit, 
EFP has a unit jump. As discussed in Section 3.3, concerning the scaling limit of 
the boundary correlation function G$ = F^' 1 ', its stepwise behaviour was related 
to the position of the saddle-point solution with respect to the pole at z = 1. It 

(r s) 

is easy to verify that the same mechanism holds for F N ' for any finite value of s: 
indeed, in the r, N — > oo limit, with s and r/N = 1 — x kept fixed (and, obviously, 
s/N = y — ► 0), saddle-point equations (|5.13[) decouple into s identical saddle-point 
equations of the form (|3 . 1 5[) . The line of reasoning discussed in Section 3.3 can be 
applied again, and due also to identity (|5.10|) . it is clear that the unit jump defining 
the position of the contact point occurs now when all s saddle-point solutions are 
located at z = 1. 

Before going on to the large s situation, it is worth to recall a peculiar prop- 
erty of the saddle-point solutions for generic Penner models (see [46] and [47] for 
further details). In these models the logarithmic wells in the potential can behave 
as condensation germs for the saddle-point solutions. It can be shown that conden- 
sation can occur only for charges less than or equal to 1, the value of the charge 
corresponding to the fraction of condensed solutions. Within the electrostatic anal- 
ogy, this corresponds to the capability for the logarithmic potential well created by 
external charge Q to screen exactly a fraction Q of the s particles of charge 1/s. 

In our case, the only possibility for condensation of solutions of the saddle-point 
equations is given by the charge at z = 1 in the triple Penner potential, since the 
charge at z = —1/t 2 is always repulsive, while the one at z = is larger than 1 in the 
region where -Tnw lies. Note further that in our case the charge at z = 1 is exactly 
1, allowing for total condensation of solutions of the saddle-point equations. This 
consideration, and the crucial identity (|5.10|) , strongly suggests that the mechanism 

(r s) 

producing the stepwise behaviour of ' in the scaling limit is still the same as in 
the case of finite s, and thus that the curve -Tnw occurs in correspondence to the 
condensation of (almost) all solutions of saddle-point equations (|5 . 13[) at the point 
z = l 

In the framework of the Green function approach, our claim can be rephrased 
as follows: the curve Tnw can be derived from the condition that x and y should 
be such that in the scaling limit Green function (15.141) reduces to 



Indeed, in [37) it is shown that this requirement translates into the following con- 
dition: 



G(z) 



1 



(5.15) 



z-1 



(1 + t 2 )(x + y - l) 2 + i 2 (l + t 2 )(x - yf = t 2 



(5.16) 
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for x,y € [0, «]. Recalling how the remaining portions of the arctic curve can 
be reconstructed (see discussion in Section 2.3), we also find that equation (|5. 16|) 
describes the whole curve A, i.e., it is valid for all x,y £ [0,1]. At A = tt/2, or 
t = 1, it turns into (a; — 1/2) 2 + (y — 1/2) 2 = 1/4, which is the arctic circle of pQ. 



5.3. Condensation of roots and arctic curve. As just discussed, the arctic 
curve 7~nw occurs in correspondence to the situation where almost all roots of 
saddle-point equations (|5.13[) condense at point 2 = 1. It turns out that this 
observation can be used to formulate a simpler approach, which does not involve 
the Green function, and can be extended to more general situations (beyond the 
case A = 0). 

The origin of the observed correspondence between the arctic curve Inw and 
the condensation of almost all roots of saddle-point equations relies on the two 
following important ingredients. On one hand, by construction, EFP in the scaling 
limit has a stepwise behaviour, the unit jump occurring in correspondence to /nw- 
On the other hand, generalized Penner models (whose saddle-point equations share 
some essential features with the ones under consideration) allow for partial or total 
condensation of roots. The role of the unit charge potential in a generalized Penner 
model with total condensation is played in our case by the s poles at z = 1 in 
representation (|5.4|) . which all turn out to be poles of order s. 

Having in mind these essential features, the correspondence between the arctic 
curve and the condensation of roots can be seen as the consequence of the following 
properties of multiple integral representation for EFP: 

• in each integration variable z\, . . . , z s the integrand has a pole at z = 1, 
and the cumulative residue over all these variables at this pole is equal to 

i; 

• the pole at z = 1 for each integration variable Z\, . . . , z s is of order s. 

The first property means that the scaling limit of EFP is governed by the position 
of the roots of saddle-point equations with respect to the pole at z = 1, and, 
furthermore, identity (|5. 10|) implies that deforming the integration contours trough 
this pole produces exactly the expected stepwise limiting expression for EFP. The 
second property means that point z = 1 is the point where the condensation of 
almost all roots of the corresponding saddle-point equations may occur. 

The properties listed above, together with the expected limiting expression 
(14. 4[) . implies that the arctic curve can be found through the condition of conden- 
sation of almost all roots of the saddle-point equations at z = 1. Let n c and n u 
denote the number of condensed and uncondensed roots, respectively, n c + n u = s. 
Condensation of 'almost all' roots, simply means that 

— ~0, s^oo. (5.17) 
s s 

Note that condensation of almost all roots at z = 1 evidently implies formula 
(|5.15p and also that RHS of (|5.13|) reduces to 2/(zj — 1), at leading order for large 
s. Hence, the system of saddle-point equations simplifies to a system of n u identical, 
decoupled, equations. This is in fact one single equation determining the position 
of the n u uncondensed roots. This equation can be called 'reduced saddle-point 
equation'. It turns out that the solutions of the reduced saddle-point equation, due 
to condensation itself, must obey certain property. 
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According to the standard picture, described in [47], the solutions of the saddle- 
point equations build up cuts in the large s limit, and these cuts move in the complex 
plane, as the parameters of the model (in our case x and y) are varied. Specifically, 
in our case, when (x,y) G D (see figure [TJ, there is a cut, with complex conjugate 
end-points (let them be z a and Zb), which lies to the left of the singularity at z = 1, 
and intersects the real axis on the segment (0, 1). When we move {x,y) closer to 
the frozen region, Fnw; the end-points of the cut, while still complex conjugated, 
get closer to the real axis, with real part greater than 1. When (x,y) reaches 
the arctic curve, the end-points of the cut join at some value w on the real axis, 
w > 1, with the cut entangled around the singularity at z = 1 . This corresponds 
to the total condensation of solutions of the saddle-point equations. When (x,y) 
enters the frozen region Fnw, the end-points separates, still sticking to the real 
axis, generating now a cut to the right of the singularity at z = 1. 

Concentrating on the case of total condensation, the cut goes now from z a = w 
just above the real axis to the singularity at z = 1, surrounds it and comes back 
to Zb = w, running just below the real axis. The two portions of the cut which lie 
one in front of the other just above and below the real axis somehow compensate 
(the corresponding densities of roots cancel each other), and one is left with a pole 
at z = 1, describing the total condensation of roots. Hence, in this case, among 
the uncondensed roots, it is necessary to have a pair of coinciding roots, with value 
w, corresponding in fact to the end-points of the would-be cut, now collapsed to 
the pole at z — 1. Such pair of uncondensed root, can in principle lie anywhere on 
the portion of the real axis to the right of the condensation pole at z = 1. Their 
position, i.e. the value of w £ [l,oo), obviously depends on the position of point 
(x, y) on -Tnw- Thus the value of w naturally parameterizes the curve -Tnw from 
the top contact point (k, 0) to the left one (0, k). 

Summarizing, we arrive at the following alternative recipe for the derivation 
of the arctic curve: on the basis of the two properties of the multiple integral 
representation for EFP listed above, we assume that they are determined by the 
condition of condensation of roots. Then, given the system of saddle-point equations 
we impose condensation of 'almost all' its solutions, obtaining a reduced saddle- 
point equation, determining the position of the uncondensed roots. We then require 
the existence among them of two coinciding roots with value w € [l,oo), which 
parameterizes the points on the arctic curve. 

To illustrate and verify this recipe, we come back to saddle-point equation, 
(|5.18p . and implement the condensation of almost all roots. As already explained, 
the last term of LHS of (|5.13[) in the large s limit reduces to 2/(zj — 1). As a result, 
this term combines with the first term and we end with the following reduced 
saddle-point equation 

^-^ + ^-^=0. (5.18) 

2—1 Z 1 + t z Z 

This equation determines the position of those roots which survive in the complex 
plane and do not condense at z = 1. 

We now require the existence among them of two coinciding roots. Noting that 
the numerator of (|5.18|) is of second order in z , we can simply require the vanishing 
of its discriminant, thus obtaining (I5.16|) . Equivalently, denoting LHS of (|5.18|) by 
F(z), we can require it to be of the form F(z) = (2 — w) 2 F(z) 1 with F(z) regular 
in the vicinity of w. This translates into a system of two equations F(w) = 0, 
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F'(w) = 0, linear in its unknowns x and y, with the solution 

t 2 w 2 t 2 (w 2 -2w + l) r , , cin , 

V= ui\ TTTIi 2 I iT ' tue[l,oo). (5.19) 



i 2 w 2 + l' a (i 2 + l)(< 2 w 2 + 1) : 

It can be easily seen that, as w varies over the interval [1, oo), the point (x, y) above 
indeed describe the curve Fnw, between the two contact points (k, 0) and (0, k). 
Elimination of parameter w from expressions (|5.19|) leads again to equation (|5.16p 
for the arctic ellipse. 



6. Arctic curves in the disordered regime 

6.1. Condensation hypothesis. In the previous Section we have shown that 
the assumption of the correspondence between condensation of roots and the arctic 
curve allows one to derive this curve in a parametric form, through the requirement 
of a pair of coinciding roots for the reduced saddle-point equation. For the case 
considered there, the correspondence between condensation of roots and arctic curve 
can be verified by direct inspection of the form of the Green function, which indeed 
reduces to expression (|5.15[) for (x,y) € -Tnw- 

It turns out that in the general case, when the value of A is arbitrary, the 
method explained in Section 5.3 remains applicable! This observation is based on 
the fact that the two properties of the integrand of multiple integral representation 
for EFP listed in Section 5.3 appear to be totally independent of the value of A. 
This means that we can again assume that the arctic curve arises in correspondence 
to condensation of almost all roots of the saddle-point equations. It is to be em- 
phasized, that in the general case however such correspondence cannot be verified, 
not even a posteriori (at least by the methods at our disposal) , hence we call it here 
'condensation hypothesis'. 

Running a few steps forward, it worth to mention here that the only limitation 
which will restrict the final result to be valid only in the disordered regime (| A| < 1) 
comes from formula (|3 . 1 1[) which will play an important role below in the analysis 
of the reduced saddle-point equation. For this reason we somehow restrict ourselves 
here to the disordered regime although the discussion is in fact rather general. 

As just said, the 'condensation hypothesis' is strongly supported by the form 
of multiple integral representation (|4.10[) for EFP in the general case, satisfying 
the two crucial properties listed in Section 5.3. To show that the first property is 
indeed fulfilled, we have to consider, for r, s = 1, . . . , N, the integral 

(r , s)= (-iy(s+w Zs r r * [(t 2 - 2M) Z] + 1]- 1 

N S !(27a)«a«(«-Dc- J c - J c - f = \ zfcj ~ l) s 
s 1 

x n t 2 _ 2Ai + 1 n to-*) 8 

j,k=l J J l^j<k^s 

x h N>s (z-i, . . . , z s )h StS (ui, . . . , u s ) dzi ■ ■ ■ dz s . (6.1) 

This integral differs from (|4.10p only in the integration contour (recall that C]~ 
denotes the closed negative-oriented contour in the complex plane, enclosing point 
2 = 1, and no other singularity of the integrand). To show that quantity (|6.1| is 
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identically equal to one, let us resort to another, essentially equivalent, representa- 
tion 



Ar lS ) ._ (-l) s / / A [ft 2 - 2At) Zj + 
V • (2*i)'& ^ fr [ \ zriz, - iy i • > 



n 



i 2 *,-^ - 2Aiz, + 1 



— - h N , s (zi, . . . ,z a ) dzi • • -dz s . (6.2) 



It is to be emphasized that equivalence of representations (|6.ip and (|6.2p simply 
follows from the fact that the integrand of (|4.10[) is the symmetrized version of the 
integrand of (|4.8p (see [3T] for details, and also Appendix C of [H] where a key 
identity has been proven). Performing integration in representation (|6.2p in the 
variable z s , and recalling property (|4.6p . the identity 



Ip s) = l, (r,s = l,...,N) (6.3) 

follows immediately by induction. In this way we observe that the first property, 
which states that the cumulative residue of the integrand for EFP over all variables 
z%,...,z s at the pole z = 1 is equal to 1, is indeed fulfilled for arbitrary value of A. 

The second property, that the pole at z = 1 for every integration variable 
zi,...,z a is of order s, follows by simple inspection of representation (|4. 10|) . This 
means that we can again expect condensation of almost all roots of the saddle-point 
equations at z = 1. Note that the other poles in the integrand of representation 
(I4.10p . being poles of small order even in the scaling limit, can only give sublcading 
contributions, with respect to the stepwise behaviour generated by condensation 
at z = 1. It is to be emphasized that although in the general situation we have 
no connection with the matrix model formulation any more, we can nevertheless 
represent the integrand in the exponential form with an 'action' which will con- 
tain, among many others terms, the term sX^=i bi(Zj — 1). Assuming validity of 
the condensation hypothesis we thus assume that this term still dominates when 
the parameters (x and y) are near the arctic curve and hence other terms have 
no relevance for the mechanism of condensation (although these other terms can 
contribute to the reduced saddle-point equation and therefore determine the arctic 
curve). 

Thus we have just seen that the two main properties of the multiple integral 
representation (|4. 10|) for EFP are fulfilled, strongly supporting the condensation 
hypothesis. Hence, we may apply the procedure explained in Section 5.3. This 
will be done in the remaining part of this section. Namely, we will derive the 
corresponding reduced saddle-point equation, and implement the condition of a 
pair of coinciding roots, to obtain the arctic curve in a parametric form. 



6.2. 'Reduced' saddle-point equation. We thus assume condensation of 
solutions of saddle-point equations (14.12|) . and derive the corresponding reduced 
saddle-point equation. We start with setting n c of the s variables z^, k = 1, . . . , s, 
to the value 1. We are left with a system of n u equations in the n u uncondensed 
variables, let them be Zj, j = 1, . . . , n u . In what follows, the fact that n u /s vanishes 
in the scaling limit, see (|5.17|) . plays a crucial role. 

In saddle-point equation (|4. 1 2(1 . let us consider the last term in the sum. As 
explained in Section 5.3, condensation of almost all roots reduces it to 2s/ {z,j — 1), 
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at leading order for large s. Similarly, for the remaining two terms in this sum, for 
large s, we have, respectively, 



t 2 z k - 2At (t 2 - 2At) 



fe= 



and 



J t 2 z 3 z k - 2Atzj + 1 (t 2 - 2At)zj + 1 



E t 2 z jZk - 2Aiz fe + 1 ~^ Vzj - 2A< + 1 + (6 ' 5) 



•o(a), (6.4) 



k=l 



Further, due to property (|4.6[) function /ij\r )S (zi, . . . , £ s ), appearing in the term 
after the sum in equation (|4.12[) , simplifies to function h N,n u {%i > • • • > z n u ) , which, 
in turn, for N, s — > oo, and n u /N ~ 0, can be evaluated for large s directly from 
its definition (|4.5[) . In this way we obtain 



\nh NjS (zi, ...,z s ) — * lnh N ,n u ( z i) ■ 



E 



lnhpf(zj) + o(s). (6-6) 



Similarly, recalling that uu — ► as — > 1, sec (|4.1ip . and using property (|4.7[) . we 
find that function h SyS {u\, . . . ,u s ), appearing in the last term of equation (|4.12[) . 
simplifies, modulo an unessential factor, to function h nutnu (ux, . . . , u nu ). Recalling 
that hn ut n u (ui, . . . ,u nu ) is a polynomial of order n u in each of its variables, we 
obtain that its logarithm for large s is estimated as o(s). Thus we have 

s 

lnh s>s (ui,...,u s ) — >liih nu , nu (ui,...,u nu )+ ^ \nhj(0) 

= C lS 2 + C 2 s + o(s). (6.7) 

Here C\ and C2 are some quantities which do not depend on Uj (j = 1, . . . ,n u ). 
After differentiating, we are left with a term estimated as o(s), which is to be 
neglected, at the leading order in the large s limit. 

As a result, we obtain that saddle-point equations (|4. 1 2[) simplify to a set of 
n u identical and decoupled equations, 

F(z j )=0 > (j = l,...,n u ). (6.8) 

Here function F(z) is given by 

F[z) := _v — izf - _ y {2 A — + ( iim ^M^y (69) 

w z-l z t 2 z-2At + l \n-,oo N J ' y ' 



where the scaling variables x and y are defined in 

To write the explicit form of the equation determining the location of uncon- 
densed saddle-point solutions, we need to know the last term in the expression for 
F(z). In what follows we shall need to consider the case of z real and positive. For 
these values we can use formula (13. for large N limit of function Hn(z)- From 
this formula, for the last term in (16.91) we obtain 



lim hi few (7(0) =]n f sina(A-77)sinagsin(g + A-77) '. (j . . |() . 

a^oo TV \asin(A — rj) sin£smct!(£ + A — rj)J 
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Taking into account that function 7(£), see (|3.12j) . is a monotonously growing func- 
tion, from — oo to +00, over the interval £ G ( — A — 77, 7r — A — 77), which is the interval 
of periodicity of 7(£) (note that in fact 7(£) is a rational function of cot£), we can 
switch from function F(z) to function /(£), defined by 

sin(A-,)siny + A + r,) 
wv " sin(A + 77) sin 2?7 J v ; 

Here the prefactor in fact is equal to l/7'(£). Direct calculation from (|6.8[) lead us 
to the following neat expression for this function 

f(0 = x^ + X)+y <p(t + 77) - (6.12) 

where function </?(£) is given by ()2.14j) and function reads 

:= cot £ - cot(£ + A + 77) - a cot(af) + a cot a(£ + A — 77). (6.13) 

Noting properties ^(^ — A — 77 — £) = and <^(7r — £) = <^(£), it easy to see that 
function /(£) obeys the symmetry under the substitution £h^7t — A — 77 — £ and 
simultaneous interchange of the coordinates, x «-> 77. This property of function /(£) 
will lead to the diagonal-reflection symmetry of the arctic curve. 

6.3. The arctic curve. To obtain the arctic curve from the reduced saddle- 
point equation, we will follow the recipe discussed in detail in Section 5.3. We thus 
have to require that function F(z) has a double zero, which, moreover, must lie 
on the interval [l,oo). Denoting the value of this zero by w, we therefore require 
that F(z) — (z — w) 2 F(z), with F(z) regular in the vicinity of w. This condition 
is equivalent to the system of two equations 

F(w) = 0, F'(w) = 0. (6.14) 

These are equations for unknown x and y, which are functions of w. For each value 
of w G [1, 00), values of x and y correspond to a point of the arctic curve, i.e., the 
solution of system (|6.14p is just the arctic curve in a parametric form, with w being 
the parameter. 

Taking into account the observation above that instead of function F{z) we 
can use function /(£), introduced by formula (|6.11j) . we immediately obtain that 
system of equations (|6.14|) is equivalent to the following system of equations 

/(C) = 0, /'(C) = 0, (6.15) 

where C is the new parameter which parameterizes the arctic curve, w — 7(C)- 
The range of values of the original parameter w, varying in the interval [l,oo), 
corresponds C € [0, n — A — 77]. 

We now solve the linear system of the two equations above, and arrive at the 
following result: 

y>'(C + vM0-v(( + vW(0 

X <p(C + AV(C + 77)-^(C + 77y(C + A)' 

^(C + A)r(C)-^(C + AMC) 1 ■ ' 

y ¥>(C + AMC + i7)-¥>(C + »?MC + A)' 
These two equations constitute the parametric form of the top left portion, InWi 
of the arctic curve, as the parameter C varies in the interval [0, tt — A — 77] . The value 
C = corresponds to the contact point of the curve Jnw with the x-axis, and as 
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£ increases, the whole curve Inw is monotonously constructed, up to the contact 
point with the y-axis, at ( = tt — A — r\. 

Using the properties of the functions involved here, we can write the result in 
the form 

x = X(0, y = Y(0, Ce[0 )7 r-A-77], (6.17) 
where functions X(Q = A(£; A, 77) and Y(() = Y((; A, 77) are simply related by 

X(0=Y(7r-X-T,-0- (6-18) 

Direct calculation gives 

sin 2 £sin 2 (£ + 277) sin(£ + A — 77) sin(£ + A + rf) 



no 



2-q sin(A — 77) [ sin(£ + A — rf) sin ( + sin(( + A + 77) sin(£ + 2?7)] 

sin(A — 77) sin(A + 7/) 
sin 2 £ sin(C + A + rf) sin(£ + A — if) 

sin(2£ + 2A) asina(A — 77) 



sin(C + A — 77) sin(£ + A + 77) sin aC, sin a(£ + A — rf) 

a 2 sina(2C + A-77)sina(A-?7)l 

• o > ■ o /a ^ \ (■ (6.19) 

sm^ac; sin^a(£ + A — 77) J 

Formulae (|6.17p - (|6.19|) represent the main result of the present paper. 

As a comment to the result, it is worth to mention that in all examples discussed 
in the literature to date, the arctic curves (or frozen boundaries of limit shapes) 
always turn out to be algebraic curves. From the form of function Y(Cf) above it is 
straightforward to conclude that as far as the parameter a = n/ (ir—2r]) is irrational, 
the arctic curve of the domain-wall six- vertex model is a non-algebraic curve. This 
is the situation for generic weights of the six-vertex model in the disordered regime. 
On the other hand, the arctic curve is algebraic when a is a rational number, 
since then both cot£ and cota£, rationally entering in (|6.19[) . can be expressed as 
rational functions of a suitably chosen parameter. In other words, the arctic curve 
in the domain-wall six-vertex model is rational only when the weights correspond 
to the so-called root-of-unity cases. 

7. Particular cases and combinatorial applications 

7.1. Alternating sign matrices. The arctic curve, whose general expression 
has been given in the previous Section, for arbitrary weights of the six-vertex model 
with domain wall boundary conditions in its disordered regime, is worth to be 
investigated in more detail when some of the parameters of the model are specialized 
to certain values. In this section we consider a few particular cases, some of which 
have a natural interpretation in the context of the limit shape of large alternating 
sign matrices. Indeed, our result for the arctic curve of the six-vertex model allows 
us to obtain the limit shapes of alternating sign matrices within g-enumeration 
schemes (for the disordered regime < q ^ 4, see below). 

We recall that an alternating sign matrix is a matrix which has only l's, O's and 
— l's in its entries, obeying, moreover, the rule that in each row and each column of 
the matrix all nonzero entries alternate in sign, and the first and the last nonzero 
entries are l's. In q-enumeration, denoted here as Ajv(g), each matrix has weight q k 
where k is the number of —l's in its entries. For certain values of q, namely q = 1, 2, 
and 3, the numbers of g-enumerated alternating sign matrices, ^4jv(<z), are known 
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to be given by some factorized explicit expressions. These results, stated first as 
conjectures, remained challenging for their proofs for a long time, see [I2ll32l l33] 
and references therein; the story and backgrounds of the problem can also be found 
in book [34] . 

A possible approach to alternating sign matrices exploits their close relationship 
with the domain-wall six-vertex model. Indeed, in |12j it was noticed that there is 
a one-to-one correspondence between N-by-N alternating sign matrices and con- 
figurations of the domain-wall six-vertex model on the N-by-N lattice. Due to 
this correspondence, An(q) is equal, modulo a simple factor, to Zjq, the partition 
function of the domain-wall six-vertex model. The weights of this model have to 
satisfy the relations a = b and, since c-weights comes in pairs, c = a^/q. The first 
relation is fulfilled when A = 7r/2 in (|2.6[) , and the second one implies that q and 
A (recall that A = cos 2rj) are related by 

<7 = 2-2A. (7.1) 

The precise relation between g-enumeration of N-by-N alternating sign matrices, 
A^(q), and the partition function Zn (at A = 7r/2) is 

Z N = a N2 q N l 2 A N {q), (7.2) 

see, e.g., [12] for a discussion. 

In the context of alternating sign matrices, the arctic curve of the domain- 
wall six- vertex model describes what is usually referred to as the 'limit shape'. 
Indeed, in their corner regions, alternating sign matrices mostly contain O's, while 
in the interior they have many nonzero entries; as the size of matrices increases, 
the probability of finding l's and — l's in the corner regions vanishes, while in the 
central region such probability remains finite 35J. As a result, in the scaling limit 
(i.e., when large matrices are scaled to a unit square) the arctic curve of the domain- 
wall six-vertex model precisely describes the shape of this central region (which is 
nothing but the region D, see Section 4). 

A derivation of the limit shape of large q-cnumcratcd alternating sign matrices 
at q = 1 and q = 3 (in addition to the well-known case of q = 2) was given in our 
previous paper [36]. There we used the observation of [41] that at r\ = tt/6, it/4:, 7r/3 
(with A = 7r/2), corresponding to the cases of q = 1,2,3, respectively, the Hankel 
determinant standing in formula (|2.13p for Zjq turns out to be related with cer- 
tain classical orthogonal polynomials (this explains, additionally, the 'roundness' 
of Ajv(<?) at q = 1,2,3). At these values of parameters the function /ijv(z) turns 
out to be expressible in terms of hypergeometric series and its large N asymptotic 
behaviour can be found directly, thus allowing us to derive the arctic curve on the 
basis of the 'condensation hypothesis'. Below we discuss several particular cases, 
and, in particular, how the results given in [36] follow from formulae (|6.17p - (|6.19l) . 
describing the arctic curve. 

7.2. Particular cases. Let us consider the arctic curve for some particular 
values of weights of the disordered regime. As mentioned above, if the parameter 
A is specialized to n/2, then the arctic curve also describes the limit shape of large 
alternating sign matrices, within the corresponding q-enumeration scheme. We 
consider here the cases which correspond to q = 1,2,3, and q = 4 where this last 
case is obtained as a limiting A — > — 1. We also consider another limiting 

case, where q vanishes. Interestingly, in this case the limit shape tends to certain 
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nontrivial limiting curve, rather than becoming somehow degenerate or trivial; such 
limiting curve is discussed in more details below. 

7.2.1. The case A = 0. We start with the case A = 0, or 77 = 7r/4. In this case, 
of course, we expect to reproduce the arctic ellipse discussed in detail in Section 5. 
Indeed, setting r\ — ir/4 in formula (]6 . 19[) . we obtain 

Y(Q = \{1 -cos 2C), Ce[0,2f-A]. (7.3) 

We also have = |(1 + sin(2£ + 2A)), and we can easily eliminate parameter 

£ in equations (|6.17p . As a result, we arrive again to equation (|5.16[) . where t := 
tan(A - 7r/4). 

7.2.2. The case A = |. The case A = 1/2, or 77 = tt/6, is interesting since at 
A = 7r/2 the model is equivalent to the enumeration of alternating sign matrices 
(with q — 1). Specifying 77 to the value 7r/6 and setting A = tt/2, we find that 
expression (|6.19p simplifies to 

y(C) = i-cosC, Ce[o,f]. (7.4) 

Correspondingly, we also have X{C,) = 1 — cos(7r/3 — £) and, eliminating parameter 
C, it can be found that the curve Jnw (the top-left portion of the arctic curve) is 
described by the equation 

(2x - l) 2 + (2y - l) 2 - 4xy = 1 , z,ye[0,±]. (7.5) 

This curve describes the limit shape of large alternating signs matrices [36 . Interest- 
ingly, in comparison with the arctic circle, given by equation (2y— l) 2 + (2x— l) 2 = 1, 
it has just single additional term —Axy in LHS. The property of the curve Jnw to 
be given by a quadratic equation in the rj = n/6 case holds only at A = n/2. Indeed, 
just specifying rj — n/6 but keeping A generic, one finds from (|6.19p that function 
y(C) m this case is rather bulky; the curve Jnw turns out to be described by a 
tenth order equation. 

7.2.3. The case A — —i. The case of A = —1/2, or rj — 7r/3, will be treated 
here only at A = 7r/2. At this choice of parameters the model is equivalent to q- 
enumeration of alternating sign matrices with q — 3, which is a well-known example 
of tractable enumeration (together with the cases q = 1,2). Specifying in (]6 . 19[) 
rj = 7r/3 and A = tt/2, we obtain 

"sin(^ — C) tanC" 



Y(C) = 4 



1 + 2cos2C 



l]H J2<, ' S ^" - V "^ C€[0,|]. (7.6) 



6 — 3 cos 2C — y/3 sin 2^ ' 

We also have X(Q — Y(ir/6— £). For further analysis, it is convenient to rewrite the 
parametric formulae for the curve in terms of rational functions of a suitably chosen 
parameter. For example, choosing the parameter w = sin(7r/6 + C)/ sin(7r/6 — £), 
one obtains formulae for the curve given in . Further, excluding the parameter 
w one can find the corresponding algebraic equation for the curve, which appears 
to be of the sixth order (see [35], equation (14)). 

7.2.4. The case A = — 1. This case can be obtained as the limit from the 
disordered regime; in fact, the case of A = —1 is often regarded as belonging to 
the disordered regime since the model is still disordered (at A = — 1 the model 
undergoes an infinite order phase transition). Denoting by v the rapidity variable 
of the model at A = —1, the standard parameterization of the weights of this model 
reads 

a = l + v, b = l-v, c = 2 (-Ku<l). (7.7) 
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To fit this parametrization, one can perform the scaling in the parameters by taking 
the limit 6 — ► 0, upon setting (see equation (|2.6| ) 

A = \ - vS, n = l~6. (7.8) 

Correspondingly, we also set C — p8 an d then equations (|6.17[) for the arctic curve 
will read x = X (p) and y = Y(p) , where p is the new parameter which parameterizes 
the curve, p £ [0, 1 + v], and 

X(p) = limX(pS), Y(p) = lim Y(p6), (7.9) 

where 5 also enters A and r\ as given by (|7.8p . Functions X(p) and Y(p) also depend 
on v as a parameter and, due to (|6.18p . they satisfy X(p) =Y(l + v — p). Explicitly, 
function Y(p) reads 

(2 - n) 2 { cos 5u 



4(l-w)[(l + i;)(l-p)+p 2 ] | v ; sinfpcosf(p-w) 

9 or , son COS-^W COS^(2p — f) ] 

-ttV 1-fc-t) 2 h V±t T > (7.10) 

^ L ^ M 4sin 2 fp cos 2 f (P~ v ) J 

where p € [0, 1+v]. Clearly, in this case, contrarily to the three examples considered 
above, the arctic curve is not an algebraic one. In the case of v = 0, that is, when the 
weights satisfy a = 6, the arctic curve describes the limit shape of large alternating 
sign matrices within the ^-enumeration at q = 4. 

7.2.5. The case A — * 1. This case corresponds to r\ — * 0. For arbitrary A and 
small r) formula (|6.19[) reads 

K(0 = ^-(2C-Bin2C) + O(»j), Ce [O.tt-A-t/]. (7.11) 

Correspondingly, we also have X(£) = Y"(7T — A — rj). The subsequent analysis 
depends on whether or not, as rj vanishes, the parameter A is scaled accordingly. 

Namely, the first possibility is to set A = rjv or A = tt — rjv where v is a new 
rapidity variable (v ^ 1); this choice corresponds to approaching either of the two 
branches, a > b or a < b, respectively, of the model at A = 1. The weights of 
the model at A = 1 are parameterized as a = v ± 1, b = v =p 1, and c = 2, where 
v ^ 1. The two choices of the signs corresponds to the two branches of the model 
at A = 1. At a > b, we find, after eliminating the parameter Q, that the curve Inw 
is just the straight line: x + y = 1, where x, y € [0, 1]. At a < 6 the curve Inw is 
just single point: x = y = 0. All this is in agreement with the fact that at A = 1 
the model is not disordered anymore; the region D, see figured! degenerates into 
the straight line x + y = 1 (if a > b) or x — y (if a < b) . 

The second possibility, which appears to be also the most interesting, is to keep 
A fixed as r\ vanishes. In the phase diagram of the six-vertex model in the a/c—b/c 
plane this corresponds to taking the limit into the deep infinity of the disordered 
region, rather than approaching either of the two branches of the A = 1 model. In 
particular, setting X — tt/2 and neglecting small rj corrections, one can easily find 
that the curve -Tnw is given by the equation 

x + y = coS7r(a;- y), x,j/e[0, |]. (7.12) 
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Figure 2 . The limit shapes of the large alternating sign matrices 
within various ^-enumeration schemes. As q vanishes, the limit 
shape tends to a limiting curve (dashed line), given by equation 

(Eli. 

This equation has an interesting meaning in the context of alternating sign matrices: 
it is the limiting curve describing the limit shape as q tends to zero. 

7.3. Discussion. A natural question concerns the qualitative behaviour of 
arctic curve (|6.17p - (|6.19p as one varies parameters of the model. In discussing 
some properties of the arctic curve it is useful to resort to numerical plotting. For 
example, considering the case of generic A, we just mention here that as A varies, the 
curve is deformed along one of the two diagonals according to the sign of A — tt/2. 
An example of the whole arctic curve A in a non-symmetric case is shown in figure 
[T] Namely, this figure shows the arctic curve at r\ = 7r/6 and A = 57r/12. 

Specializing to the case of A = tt/2 we can also discuss the arctic curve in 
application to the limit shape of large g-enumerated alternating sign matrices, with 
q and rj related by (|7.1[) . The relevance of formulae (|6.17p ~ (|6.19p in the context of 
alternating sign matrices is that they allow one to study the variation of the limit 
shape as q varies over the interval (0,4]. In addition to the cases of q = 1,2, 3 
considered in [3S] , in figure we plot also the limiting cases of q = 4 and q — > 
(the latter shown by a dashed line). Note that, as expected, the disordered region 
(the area enclosed by the limit shape) is largest for q — > 0, and slowly decreases as 
q increases over the considered interval. 

Comparing numerical plots of the arctic curve at various values of the pa- 
rameters of the model we find that our results are completely compatible with all 
available numerical data [T51[T5], which are however affected by large uncertainties, 
due to the huge technical difficulties in this kind of computer simulations. The 
most refined computer simulations available at the moment has been performed by 
Ben Wieland for the especially relevant case of alternating sign matrices at q = 1 , 
corresponding to the value A = 1/2. Pictures comparing these numerical data with 
the corresponding arctic curve, given by equation (|7.5p . are available [4"8] . 
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Coming back to the qualitative behaviour of the arctic curve, we would like 
to focus again on the limiting curve, which, in the context of limit shapes of large 
alternating sign matrices, is referred above as the 'q — > 0' curve. Note that such 
a curve exists for any fixed value of A, as r\ — > (see discussion in Section 7.2.5). 
Concerning the emergence of a non-trivial limiting curve, in hindsight, it is clear 
that this is ascribable to the fact that the two limits, N — > oo, and q — > 0, do not 
commute. 

8. Conclusions 

In the present paper we have derived the arctic curve for the six-vertex model 
with domain wall boundary conditions, in its disordered regime. The derivation 
is essentially based on the exact expression in terms of a multiple integral for the 
emptiness formation probability, a correlation function which discriminates order 
and disorder. We have observed that in the scaling limit the arctic curve can 
be obtained as the condition that almost all solutions of the system of coupled 
saddle-point equations for the multiple integral representation of emptiness forma- 
tion probability condense at the same, known, value. The explicit expression for the 
curve in parametric form shows that in general it is a non-algebraic curve; for spe- 
cial choices of weights of the model (the so-called root-of-unity cases) it simplifies 
to algebraic curves. We have also discussed combinatorial applications of the result 
to the problem of limit shapes of large alternating sign matrices within various 
enumeration schemes. In particular, we find that the limit shape of q-enumerated 
alternating sign matrices has a non-trivial limit as q tends to zero; furthermore, 
such limiting curve is described by a very simple equation. 

Having established the expression of the arctic curve, a natural question con- 
cerns its fluctuations, which, in our approach, are related to the subleading correc- 
tions to the stepwise behaviour of emptiness formation probability in the scaling 
limit. In the case of domino tiling of Aztec diamond the fluctuations of the arctic 
circle are governed by the Tracy- Widom distribution and the Airy process [MlUB] . 
This results, which naturally extends to the case A = of the domain-wall six- 
vertex model |16j . appears rather natural in view of the 'random matrix model' 
derivation of the arctic circle, provided in [37]. Indeed, the fluctuations of the arc- 
tic curve are expected to be governed by the fluctuations of the first eigenvalues 
evaporating from the logarithmic well where total condensation occurs, and these 
fluctuations are in turn known from the literature on Penner models (see, e.g., |47j ) 
to be governed by the Tracy- Widom distribution. From the discussion of Sections 
5 and 6, this pictures extends rather naturally to generic values of A in the disor- 
dered regime. On the basis of universality considerations, it is thus very tempting 
to argue that fluctuations of the arctic curve are still governed by the Airy process, 
at least in the disordered regime. A direct proof of this statement would be of great 
interest. 

Another natural question concerns the extension of our results to the anti- 
ferroelectric regime, A < — 1. We recall that in this case (see Section 2.3) there 
are two different phase-separation curves, an outer one, which is the usual arctic 
curve, separating regions of ferroelectric order from an intermediate region of dis- 
order, and an inner one, separating this region of disorder from a central region of 
anti-ferroelectric order. Since the emptiness formation probability detects spatial 
transition from order to disorder, the present approach can be used to address the 
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problem of the outer phase-separation curve. Most of its ingredients arc indepen- 
dent of the value of A, and thus it can be directly applied to the anti- ferroelectric 
regime. In particular, in the special case of A — > — oo (which is technically similar 
to the free-fermion case) the arctic curve can be readily derived, reproducing the 
result of paper 20.]. For generic values of A in the anti- ferroelectric regime, the 
only open problem concerns the evaluation of the thermodynamic limit of func- 
tion htf{z). Such evaluation would provide the solution to the problem of the 
outer phase-separation curve of the domain-wall six-vertex model for the whole 
anti-ferroelectric regime. 
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Appendix A. Partition function of partially inhomogeneous model and 

functions hjf <a {z\, . . . , z s ) 

Formula (|2.12p is a special case of a more general determinant representation, 
known as Izergin-Korepin formula, which was originally derived for the model with 
inhomogeneous weights [25] . The weights are made position-dependent by attaching 
rapidity variables to each vertical and horizontal line, so that there are 2N rapidity 
variables in total instead of just one variable A. Namely, the weights of the vertex 
lying at intersection of ith vertical and kth horizontal lines (we enumerate lines 
from right to left and from top to bottom) are parameterized as 

a ik — sin(A, i - v k + 77), b ik = sin(Aj - v k - 77), c 4fc = sin2?7. (A.l) 

Correspondingly, the partition function is now a function of 2N rapidity variables 
Ai, . . . , Ajv, v\, . . . , vjj. In [25] the following representation was shown to be valid 

nffe=i sin(Ai - v k + 77) sin(A l -v k -rf) 

xdetM, (A.2) 

where M is an N-by-N matrix with the entries 

, , \ sin 2n , . . 

M ik := ^(A, - u k ) - — ' -. (A.3) 

sm(Aj - vk + T)) sin(A, —v k — r)) 

To obtain from (TA~2]) . one has to set v k — and = A (i, k — 1, . . . , N). 

Due to the singularities in the denominator of representation (|A.2[) . this has to be 
implemented as a limit, to be evaluated using l'Hopital's rule. 

In performing the limiting procedure to the homogeneous model there are some 
interesting intermediate situations. An important example is when the limit is 
performed only in one set of variables, say, v k — ► (k = 1, . . . , N) while all Aj's are 
kept generic (and different from each other). Therefore, the weights are given by 

ai = sin(Ai + 77), bi = sin(A^ — 77), Ci=sin2^. (A. 4) 
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and the partition function Zn(Xx, . . . , Xn) '■= Zn(Xi, . . . , Xn] 0, . . . , 0) reads 



Zjy(Xi, . . . , Xn) 



¥>(Ai) 

^(Ai) . 



{N-l) 



(Ai 



•r 



tp(X N ) 

v'(x N ) 

(Jv-D(A 



(A.5) 



This case can be called partially inhomogeneous model: the model is homogeneous 
along one direction, but still inhomogeneous in the other one. 

Having in mind this case as the starting point, one can consider, for any chosen 
s (s = 1, . . . , N), the case where Ai, . . . , X s are different, but A s+ i = • • • = A^ = A. 
One can still refer to such a situation as the partially inhomogeneous model. It 
turns out that in this case, the partition function is closely related to function 
h,N,s(ui, . . . , u 8 ) given by (I4.5|) : in particular, the partition function (|A.5[) is related 
to function /i/v,iv(wi, ■ ■ ■ ■, Un). 

Let us define variables £i, . . . ,£jv such that 



Ai=A + &, i = l,...,N. 
It is convenient to introduce the 'bare' partition function, 

Zjv(Ai, A2, . . . , Xn) 



N) 



(A.6) 



(A.7) 



Zn(X, A, . . . , A) 

Below we shall treat £j's as variables while A is to be regarded as a parameter of 
the homogeneous model. 

It is useful to consider first the case where all £,-'s are zeros but one, say £1, is 
kept nonzero. Denoting £ := £1 we straightforwardly have 

N 



Zjv(£,0,...,0) = 



(N-iy. 



(sin^^- 1 
where, for later use, we have denoted 



y(A) 



Dn(0 



D 



N 



(A. 



Dn(0 



tp'(X) 



(JV-2) 
{N-l) 



(A) 
(A) 



^ 2W - 3 )(A) ^-i)(A + £) 



(A.9) 



and Dn is given by (|2.13[) . Note that in the case of one nonzero inhomogeneity one 
can always assume that it is attached to the last column (since Zn(X\, . . . , Xn) is a 
symmetric function in its variables); thus we have chosen £2 = • ■ ■ = £v = above. 

To proceed further let us come back to the definition of the partition function 
as a sum over all configurations. The peculiarity of the domain wall boundary 
conditions is that they admit only one vertex of weight c in the last column; if this 
vertex is at rth position (counted from the top) then the first (r — 1) vertices are 
of weight b while the remaining (N — r) vertices are of weight a. As it is has been 
explained in Section 3, the probability of having the c-weight vertex at rth position 

(r) 

on last column is equal to H N , the correlation function which is originally defined 
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as the probability of having the c- weight vertex at rth position on the first row, see 
(|3.1|) . Due to this observation, the following relation is valid [33] , 



N 



Zn(Xi, A, . . . , A) — Zn ^2 H. 



(r) 
N 



sin(Ai 


- v)~ 


r—1 


sin(Ai ■ 


+- vY 


N-r 


sin(A 


-v). 




sin(A - 


-v) . 





(A.10) 



where Zm — Zn(X, . . . , A). Recalling that h^{z) — Y^r=i Hn' zT 1 an d denoting 

(A.ll) 



r(»\)»,r— 1 



sin(A + ?7) sin(A + £-77) 
7(?) = -r 



sin(A — rj) sin(A + £ + 77) 
we obtain that formula (|A.10[) actually reads 

, -1 JV-l 



Z JV (£,0,... ) 0) = 



sin(A + C + 7?) 



Mt(0). 



(A.12) 



sin(A + 77) 

More generally, as it has been first observed in [33] (at the special case A = tt/2) 
and proven in [31j . for all s = 1, . . . , N, one has 

"sin(A + & + i])' 



Zn(£i, •••>Cs,0, . 
where /ijv,s(ui) 



-o)=n 



i=l 



sin(A + 77) 



AT-1 



llN,s(ui, . . . ,U a ) , (A. 13) 



, u s ) is defined in (14. 5|) . and 

tti:=7(6)- (A- 14) 

The proof of (|A.13|) is based on formula (|A.5[) and some standard facts from the 
theory of orthogonal polynomials. We refer for further details and proofs to [3"T] . 

Finally we mention that formulae (|A.8[) and (|A.12[) provide a representation 
for generating function hjq{z). This representation is used in Appendix B for a 
derivation of the large N leading term of In hj^(z). 

Appendix B. Large N limit for function h^{z) 

Our aim here is to show how formula (|3.11|) for the large N limit of function 
/ijv(z) can be derived. To this aim we extend here the method of [27], where the 
leading term of In Zn have been found. It is to be stressed that we assume that the 
weights correspond to the disordered regime only (see |28j for a discussion of why 
the same approach cannot be used for the anti- ferroelectric regime). 

The approach is based on making use of some differential equations which can 
be obtained by means of the so-called Sylvester determinant identity. This identity 
relates the determinant of a matrix with a determinant of some other matrix whose 
entries are minors of the original matrix. Namely, let A be an n-by-n matrix, with 
entries Aj^- Let us consider minors of degree p (1 $J p $J n) of this matrix, 



ll,l2, ■ 

k\, k 2 ■ 



^*2, fcl 



A 



■^■i 2 ,k 2 

Ai„h„ 



i 2 ,k p 



where 1 ^ i\ < 12 < 
matrix B such that 



< i p $C 71 and 1 k\ < k% < 



■ ■ Ai k 

< k p < n. 



(B.l) 



Let us introduce 



B, 



k •= 



A 



1,2,. ..,p,p + i 
1,2,. ..,p,p+k 



i,k = 1, . 



(B.2) 
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Hence, entries of B are minors of matrix A of degree (p+1). The Sylvester identity 
reads 

—(n—p—l) 



detA 



A 



detB. 



(B.3) 



1,2,...,?' 
k l,2,...,p, 

The standard proof of this identity is based on the Gauss algorithm, see, e.g., 
monograph [25] for further details. In all examples below we put n = N + 1 and 
p = N — 1, i.e., B will be a two-by-two matrix. 

Let us consider first the partition function, Zjy, given by (|2.12|) . We start with 
mentioning that for arbitrary N the partition function Zjy = Zpj(\) satisfies 



Z N (X) = Z N (n - A), Z N {X)\ x=r] = (sin 2??) 



(B.4) 



Here, the first relation reflects the crossing symmetry of the model (since 6(A) = 
a(ir— A)), and the second one follows from the fact that if b = (or a = 0) then there 
is exactly one configuration contributing to the partition function. As N — > oo, one 
has [27] 

Z N = exp(-N 2 f + 0(N)), (B.5) 

where / is the free energy per site, the quantity of interest. To find /, let us consider 
the Hankel determinant in (I2.12| . Dm, given by (|2.13|) . As N — * oo, we have 

-i 2 



D N = 



N-l 



exp(N 2 <t) + 0(N)) , 



where <fi and / are related by 

— / = lnsin(A + 77) + lnsin(A — rj) + <p. 



(B.6) 



(B.7) 



To compute <fi, let us consider D^+i and apply the Sylvester identity, with 
n = N + 1 and p = N — 1 , which gives 



D 



1 



JV+l 



d 2 x D N 
dxD N 



d x D N 
D N 



(B. 



In writing (|B.8|) we have taken into account that the first derivative of a Hankel 
determinant changes the last row or the last column only, and the second derivative 
changes both the last row and the last column only. Relation (|B.8|) can be rewritten 

as 



D 2 

N 



'N • 



(B.9) 



From (|B.6[) it follows that, as N — > 00, the leading terms in each side of (|B.9[) are 
of order N 2 ; marching terms gives: 



To solve this equation it is convenient to define W := exp(— 
at 

(d x W) 2 - Wd\W = 1. 
Obviously, the general solution of this equation has the form 

sin(ai A + a-i) 



W = 



Oil 



(B.10) 
so that we arrive 

(B.ll) 
(B.12) 
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where a% and a 2 are some constants. To fix these constants, let us turn to properties 
(|B.4|) . and via relation (|B.7j) we obtain that they imply W(A) = W(ir — A) and 
W"(A)| A= = 0, respectively. We find 



a 2 



Oil 



7T — 2?7 



(B.13) 



In terms of parameter a given by (|2.17p . we have a\ = a and 0,2 = —Tjot, and 
therefore the result for 6 reads 



In 



sin a (A — rf) 



(B.14) 



Substituting this result into (|B.7|) . we readily reproduce formula (|2.16p for the free 
energy in the disordered regime. 

A detailed analysis of the large N expansion for the partition function of the 
domain- wall six- vertex model can be found in [39j . 

Let us now turn to the function hpj(z). Due to formula (j A. 12|1 . the large N limit 
of function Hn(z) can be found from that of the partition function Zn(£, 0, . . . , 0). 
This quantity is given by formula (1A.8|) . where the nontrivial object is the last 
factor, the ratio of two determinants, which we denote 



D 



N 



(B.15) 



(note that all quantities here also depend on A, which is to be regarded now as a 
parameter). For any finite N we have the property (recall that Zn(^, 0, . . . , 0) |^=o = 
1) 



Sn(0 



N-l 



{N-iy: 



(B.16) 



From the interpretation of 2T/v(£, 0, . . . , 0) as the 'bare' partition function it follows 
that it can grow up at most exponentially as N increases, so that, as N — * 00, we 
have 



Sn(0 



1 



(N-l)\ 



exp{Nip(£)+o(N)). 



(B.17) 



We show now that can be found as a solution of an ordinary first order differ- 
ential equation (in the variable £). In turn, the result will imply formula (j3.1ll) . 

Let us again use Sylvester identity, with n = N + 1 and p = N — 1, applying it 
to Dn+i and to 8\Dn+i (see (|A.9jl ). that gives 



D 



N+l 



(0 



D 



N-l 



dip*® d x D N 
D N (0 D N 



D 



N-l 



d\D N (£) d 2 x D N 
D N (0 d x D N 



(B.18) 
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In terms of function (|B.15j) these two relations read: 
-Djv+i-D/v-i 



D 2 

N 



s N+1 (0 = d x s N (0 



d x S N+1 (0 + Sn+1 (0 = d e S N+1 (0 + ° Nd >^ N d x S N (Q, (B.19) 

(d x D N ) 2 - D N djD N 

+ „ f, 

The first relation here allows us to eliminate 8\Sn(0 and d\SN+i{£,) in the second 
relation, and further, using (|B.9[) and shifting N 1— >• N — 1, we obtain 

DN+ 'f N - 1 S N+1 (0 + (d x In S N (0 + Sjv-i(0 = d 6 S N (0- (B.20) 

After dividing by SW(£), both sides of relation (|B.20|) are of order N in the 
large N limit; matching terms leads to the following ordinary first order differential 
equation: 

e 2 ^ + 2d X (f> + e^ (0 = S^(0- ( B - 21 ) 
Here the function <p is given by (|B.14j) ; explicitly, the equation reads 

2 

. . ° e*^ _ 2a cot a(A - n) + e - *^ = (B.22) 

sm^a(A — rj) 

Due to (|B.16|1 we need the solution of this equation such that expV>(£) ~ 6 when 
£ — * 0. The solution reads 

asina(£ + A — ry) ' 

Recalling formulae (|A.8|1 . ()B. 15|) . and (|B. 17|) . we therefore find that, as N — > oo, 

~ , . ( wfA) sina(A — 77) sinaf \ , „ 
In Z N (g, 0, ■ ■ ■ , = N In f K > — ; /; - +o (JV . B.24 

Finally, recalling that function y(A) is given by (|2.14p and that function hj^{z) is 
related to Zjv(£, 0, • • ■ ,0) by (|A.12[) . we readily obtain expression p. lip , which is 
thus proven. 
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